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We present a self-consistent field approximation to the problem of the genetic switch composed 
of two mutually repressing/activating genes. The protein and DNA state dynamics are treated 
stochastically and on equal footing. In this approach the mean influence of the proteomic cloud 
created by one gene on the action of another is self-consistently computed. Within this approxi- 
mation a broad range of stochastic genetic switches may be solved exactly in terms of finding the 
probability distribution and its moments. A much larger class of problems, such as genetic networks 
and cascades also remain exactly solvable with this approximation. We discuss in depth certain 
specific types of basic switches, which are used by biological systems and compare their behavior to 
the expectation for a deterministic switch. 



Introduction 

Genetic switch systems are an elementary means of 
regulatory control present in every living organism. Their 
complexity and details differ, but the general mechanism 
of the expression of a given gene being regulated by pro- 
teins, is believed to be universal (Ptashne and Gann, 
2002). They are building blocks of larger regulatory el- 
ements: genetic networks and signaling cascades. The 
pathways by which these systems operate is passed on 
from generation to generation. Understanding their sta- 
bility and characteristics is therefore fundamental. A 
lot of previous work has considered a deterministic de- 
scription of genetic switches (Shea and Ackers, 1982), 
(Hasty et al., 2001). The need for a stochastic treatment 
of genetic switches due to the single copy of the DNA 
molecule and multiple protein molecules in the cell, has 
been largely recognized (Sneppen and Aurell, 2002), (Ke- 
pler and Elston, 2001). 

The most general way of accounting for non determinis- 
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tic processes is to write down the master equation for a 
given system. To define the state of the switch one must 
specify the DNA binding states of particular genes and 
the number of proteins of each type. The probability dis- 
tribution even of a single switch consisting of two genes, 
the product proteins of which act as regulator proteins 
for the system, may not be determined exactly and ap- 
proximations must be considered (Bialek, 2001), (Hasty 
et al., 2000), (Sneppen and Aurell, 2002). 
Several approaches to account for the probabilistic nature 
of chemical reactions have been undertaken, ranging from 
the Langevin description of single genes (Bialek, 2001), 
and two interacting gene switches (Hasty et al., 2000), to 
the master equation reduced to a Fokker-Planck equation 
considerations (Kepler and Elston, 2001), (Hasty et al., 
2001a). A dynamical action formulation has also been 
used (Sneppen and Aurell, 2002) to determine the life- 
times of states of the switch. A popular alternative to 
purely analytical methods, which often need to make ap- 
proximations or are limited to very simple model sys- 
tems, has been to conduct stochastic simulations of ge- 
netic switches. Two types of simulations are mostly used. 
In the first the randomness of the system is introduced by 
means of a Monte Carlo algorithm with fixed time step 



2 



(Paulsson et al. , 2000) . The second is based on the Gille- 
spie algorithm (Gillespie, 1977) to predict the probabil- 
ity of a given reaction occurring (Arkin et al., 1998). For 
single gene systems, stochastic simulations have shown 
that stochasticity in the system is responsible for the bi- 
modal probability distributions (Cook et al., 1998), ob- 
served experimentally These methods prove very useful, 
as they allow us to test the theoretical predictions on 
model systems, which might be hard to build experimen- 
tally However this approach often does not enable us to 
gain intuition or insight into the mechanisms behind the 
functioning of the system. The aim of the present work 
is to gain a better and deeper understanding of the de- 
vice physics of genetic switches. We therefore, contrary 
to many important previous discussions (McAdams and 
Arkin, 1997), (Aurcll et al., 2002), (Vilar et al., 2003) do 
not present a specific concrete biological system, but dis- 
cuss generic behavior and try to understand its sources. 
Our approximation also allows for an exact solution of a 
broad class of genetic switch systems without any further 
assumptions and with little computational effort. Hasty 
et al (Hasty et al., 2001b) present an overview of the ex- 
istent theoretical approaches. 

A popular approximation, assumes the DNA binding 
state reaches equilibrium much faster than the protein 
number state. Therefore the adiabatic approximation is 
often considered (Shea and Ackers, 1982), (Sneppen and 
Aurell, 2002), (Darling et al., 2000), allowing for a ther- 
modynamic treatment (Shea and Ackers, 1982) of the 
DNA binding state. The protein number fluctuations are 
then treated stochastically. Even before the statistical 
thermodynamics approach of Shea and Ackers (Shea and 
Ackers, 1982) using partition functions, much previous 
work assumed the DNA binding and unbinding can sim- 
ply be accounted by an equilibrium constant, since the 
relaxation timescales for equilibration of the DNA state 
are much larger than those of the protein numbers, which 
require protein synthesis and degradation to change. The 
partition function approach has also been successful at 
looking at logic gates build from switches (Buchler et al., 
2003). The adiabatic approximation is believed to hold 
true in many cases, judging by the experimental param- 
eters of biological switches (Darling et al., 2000). But 
as the experiments of, for example Becskei et al (Bccskci 
et al., 2001) show, not all switches need function in the 



adiabatic limit and the non-adiabatic limit may result in 
new phenomena. We therefore consider a wide range of 
parameter ratios in our discussion. 

In this paper we explore more fully an approxima- 
tion, previously used by Sasai and Wolynes (Sasai and 
Wolynes, 2003) for the variational treatment of the prob- 
lem, the self-consistent proteomic field (SCPF) approx- 
imation. Within this approximation one assumes the 
probability of finding the switch in a given state is a 
product of probabilities of states of individual genes. One 
can then solve the steady state master equation for the 
probability distribution of many regulatory systems ex- 
actly. We discuss the approximation and present a de- 
tailed study of different classes of genetic switches, some 
of which have never previously been considered theoret- 
ically. We consider several particular features of such 
systems, found in known switches, separately to be able 
to characterize their contributions to the behavior of the 
whole system. To be specific, starting from a symmetric 
toggle switch, we go on to compare the effects of multi- 
mer binding and of the production of proteins in bursts 
on the stability of the switch. 

The stochastic effects prove to be modest for symmet- 
ric switches without bursts, especially if the genes have 
a basal production rate. We find the deterministic and 
stochastic SCPF solutions to have similar probabilities 
of particular genes to be on and mean numbers of pro- 
teins of a given species in the cell. However in the non- 
adiabatic limit, when the unbinding rate from the DNA 
is smaller than the death rate of proteins, the probabil- 
ity distributions have two well defined peaks, unlike in 
the deterministic approximation or adiabatic limit of the 
stochastic SCPF solution. 

We also show the effect of stochasticity on the observ- 
ables becomes more apparent when proteins are produced 
in bursts. In these types of switches, the definition of 
the adiabatic limit, which was clear for the switches in 
which proteins are produced separately, is no longer sim- 
ple. Our discussion shows that the properties of genes of- 
ten analyzed in the deterministic limit, may be strongly 
influenced by stochasticity in this case. Randomness in 
a biological reaction system leads to quantitative and in 
many examples even qualitative changes from predictions 
of deterministic models. 

We also discuss the differences in the behavior of an 
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asymmetric and symmetric switch. We point to the 
mechanisms resulting in different types of bifurcations 
and show how they are influenced by noise. Within 
the SCPF approximation switches that are regulated by 
binding and unbinding of monomers, do not have regions 
of bistability. This holds true for both symmetric and 
asymmetric switches. When proteins are produced indi- 
vidually rather than in bursts, fast unbinding from the 
DNA can effectively minimize the destructive effect of 
protein number fluctuations on the stability of the DNA 
binding state. Furthermore a detailed analysis of the 
probability distributions show they have long tails and 
are far from Poissonian in both the adiabatic and non- 
adiabatic limit. We discuss the properties of the system 
in terms of clouds of proteins buffering the DNA. We 
show how fast or slow DNA binding characteristics and 
protein number fluctuations influence the stability of the 
buffering clouds leading to specific emergent behavior of 
observables. Throughout the paper a comparison is made 
between results of the exact stochastic solution with so- 
lutions of deterministic kinetic equations for the system, 
within the self-consistent proteomic field approximation. 
We establish a base of potential building blocks of more 
complicated switches and systems, such as networks and 
signaling cascades, for which an exact solution within the 
present approximation can also be obtained. A detailed 
discussion of these larger systems will be the topic of an- 
other paper. We also present limitations of the present 
style of analysis where exact solutions are not possible. 
There are two aims of this paper. The first is to discuss 
the self consistent field approximation and show that it 
has an exact solution which may be extended to a large 
class of systems. This approximation lets one deal in a 
straightforward and computationally inexpensive manner 
with the effect of random processes on genetic networks. 
The second is to discuss the many components of biologi- 
cal switches present in nature and in engineered systems, 
in the necessary stochastic framework. 

The Self-Consistent Proteomic Field Approximation 

The basic mechanism of gene transcription regulation 
in prokaryotes may be reduced to the binding and un- 
binding of regulatory proteins, repressors and activators, 



to the operator site of the DNA. If we use this simplified 
treatment, which neglects extra levels of regulation, such 
as the binding of RNA polymarase, effectively each gene 
can be described as being either in an active (on) state, 
when the repressor is unbound (activator bound), or in 
an inactive (off) state, with the repressor bound (activa- 
tor unbound) . The stochastic system of a single gene and 
its product proteins is described by the joint probability 
distribution P(n,t) — (Pi(n,t), P 2 (n 7 t)) of the number 
of product proteins in the cell n, and the DNA bind- 
ing site state: on (protein not bound)- 1, or off (protein 
bound) - 2. To conserve probability ^ n P{n,t) = 1. 
If one considers two interacting genes, the description 
in terms of a joint probability vector needs to be ex- 
tended to four states: both genes may be on, or off, or 
one of the genes may be on, the other off. If the two 
genes do not interact, as would be the case for two self 
regulatory proteins, the probability of a finding the two 
gene system in a given state, defined by both the num- 
ber of product proteins and the DNA binding site state, 
would be the the product of the states of particular genes 
Pjji (rii ,n 2 ;t) = Pj («i ; t)Pj> (n 2 ; t) . This is generally not 
true for two interacting proteins, as is the case in a ge- 
netic switch. However, as a first approximation to the 
problem, one can ignore correlations between the spaces 
of the two genes and assume the space of the switch is a 
sum of spaces of the genes that compose it. Since we are 
looking for solutions in which the symmetry of the system 
is broken and different behaviors of the on and off state 
of a gene are possible, we must allow for different prob- 
ability distribution functions for the on and off states. 
This is analogous to the unrestricted Hartree approxima- 
tion in quantum mechanics, where allowing different spa- 
tial functions for spin up and spin down states results in 
breaking of the symmetry of the bound molecular orbital 
solution to the dissociated solution of two separate hydro- 
gen atoms with opposite spin states for large internuclear 
distances. We therefore allow for multiple solutions for a 
given set of parameters. The total probability of having 
a given gene state i and rii proteins of that type is simply 
given by Pj(ni,n v ) = Pj,j>=o(ni,ni>) + P jtj/=1 (ni, n^). 
The self-consistent approximation is a crude approxima- 
tion since in the case of the genetic switch, the state 
of a given gene is determined by the number of protein 
products of the other gene. However, within this ap- 
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FIG. 1: A schematic representation of the toggle switch. Gene 

1 produces proteins of type 1 which repress gene 2 and gene 

2 produces proteins of type 2 which repress gene 1. 

proximation, one can solve the master equation for the 
probability distribution exactly without any further ap- 
proximations. This yields a powerful computational tool, 
which simultaneously gives useful insight. 



The Toggle Switch 

For clarity of exposition, we show how the problem 
may be solved exactly within the self-consistent field ap- 
proximation on a well defined system of the toggle switch. 
We then expand the method to apply to other systems. 
The elementary system we use as an example is com- 
posed of two genes, labeled 1 and 2, as presented in Fig. 
n Gene 1 produces proteins of type 1 which, act as 
regulatory proteins, say repressors, on gene 2. The prod- 
uct of gene 2, proteins of type 2, in turn repress gene 
1. In this simplified model, we assume that protein pro- 
duction occurs instantaneously upon unbinding of the re- 
pressor. For now, we assume that repressor proteins bind 
as dimers, since that is a common scenario in biological 
systems, but we do not treat dimerization kinetics ex- 
plicitly. For simplicity the coupling form between the 
genes responsible for binding will be taken to be of the 
form hin^i, where p is the order of the multimcrization 
of the repressor. This form is a small approximation to 
the more exact /ijn3_j(n3_, — l)...(ri3_i — p + 1). We 
have checked that using the simpler monomial does not 
influence the results in any regime discussed. We also 
do not account for the existence of mRNA molecules and 
the consequent time delays owing to their synthesis as in- 



termediates. The extensions of the model are discussed 
later. Within the self-consistent field approximation the 
set of master equations for the corresponding system is 
of the form: 



dPijnj) 
dt 



dP 2 (n l ) 
dt 



gi{i)[Pi{ni-l)-Pi{ni)} + 
+h[( ni + + 1) - niPxim)] + 

-h i n\_ i P 1 {n. l ) + fiP 2 (ni) 

92{i)[P2(n i -l)-P 2 {n i )} + 

+h[( ni + \)P 2 {n l + 1) - mPiim)] + 

+h i nl_ i Pi(n i ) - fiP 2 (rii) 

for n > 1 where the i — 1,2 refers to the gene label. 
P\(ni) describes the probability of gene 1 being in the 
on state and there being m protein molecules of type 
1 in the cell. The first term on the right hand side 
of the equations describes the production of proteins 
of type i with a production rate <?j(«), where j=l,2, 
depending on whether the gene is in the on or off 
state. The second term accounts for the destruction 
of proteins with rate fc^. The binding of repressor 
proteins produced by the other gene is proportional to 
the number of dimer molecules present in the system 
ri3_i with rate hi. We assume unbinding occurs with 
a constant rate fi. Binding and unbinding contributes 
to the kinetics of the DNA binding states, as described 
by the last two terms. This set is supplemented by the 
Pj(fk = 0) equations to account for boundary conditions. 



dP 1 (n i =0) 
dt 

dP 2 {n z = 0) 
dt 



-gi(i)Pi(n i = 0) + kiP 1 (n i = 1) 
-hinl^Pxim = 0) + fiP 2 {ni = 0) 
-g 2 (i)P 2 (n< = 0) + hP 2 (n % = 1) 
+h i nl_ i P 1 {n i = 0) - f.P^n, = 0) 



For convenience, let us define ^ n Pjij^i) = Cj> the prob- 
ability of finding the DNA binding site in a given state. 
One can now sum the Pj{l) equations over the number 
states of the 2-the proteins with Pi (2) + P 2 {2), and like- 
wise the Pj(l) equations. Due to the SCPF approxima- 
tion, the only term affected is the repressor binding term 
hi(n 2 ), and since J2n 2 -fi (2) + -f 2 (2) = 1, the summation 
results in£„ 2 ft 1 (n|)(P 1 (2)+P 2 (2)) = hx{Cx{2) < n\ 2 > 
+C 2 (2) < n| 2 >) = hiF(2), where < n| 2 > is the second 
moment of the number distributions of type 2 proteins 
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produced when gene 2 is in the j-th state. The equations 
of motion of the moments of the probability distribution 
are of the form: 



roth order equations one can calculate 



dC (i) < n k > 

3 m 3 ~ = 9j(i)[< fai + l) k >-< n% >]Cj(i) + 

+ki[< njiiriji - l) k > - < n^ +1 >]Cj(i) + 

+ (-iyh t F(3-i) <n k li >C 1 (i) + 

+{-i) 3+1 h < n k 2t > C 2 (i) 

The steady state equations for the moments of the distri- 
butions that follow are closed form, the nf 1 order moment 
equation of motion depends only on the lower moments 
of the i th gene and n'l_ i . 

To analyze the behavior of switches we introduce the 
following scaled parameters: the adiabaticity parameter 
Wi = fi/ki, which represents the characteristic rate of 
change of the DNA state compared to the characteristic 
rate of change in protein number, X^ q = fi/hi measures 
the tendency for proteins to be unbound from the DNA, 
Xf 1 = 5(51 (*) + 92{i))/ki the effective production rate 
and 5Xf w — \{gi(i) — g2(i))/ki distinguishes between 
the two DNA states in terms of protein dynamics. Fur- 
thermore, in the operator formalism developed for classi- 
cal diffusion by Doi (Doi, 1976) and Zeldovich' and Ov- 
chinikov (Zeldovich and Ovchinikov, 1978), the number 
operator may be written in terms of number state cre- 
ation and annihilation a operators, as n = a^a. It 
is then particularly easy to write down the equations for 
the a moments instead of the n moments. Setting the left 
hand side to zero one obtains the steady state equations: 

= - Wi [^^Cr 1 (i)-C 2 (i)] 



< 4 - 4 >= {{Xf + 6XD < a^ 1 > + 

(Y ad A ysw\ < fe-1 ^\ kCj{l) 



which depends only on a moments of lower order than 
the k th moment. This allows one to obtain the following 
form for the higher order a moments 



<a k u > 



<4i > 



(Y ad XY SW \(~I UjC^i) . k _ 1 

{Xi ~ SXi ^ i + kC 1 (i) )<<h > + 



Going back and forth between the two types of moments 
is straightforward. The n-moment equations have how- 
ever more complicated forms: 



= k[(Xf + (-lySXD < a)' 1 > - < a% >]Cj(i) + 



req 



< a k u > < 0% > C 2 {i)\ 



Using the probability conservation relation C\(i) 
C 2 {i) = 1) the zeroth order equations become: 



Ci(i) 



x 



eg 



Xt q + F(3 - i) 



C 2 (i) 



F(3 - i) 



Xl q + F(3 - i) 



(1) 



Dividing the higher order dj{i) moment equations by 
Cj(i) and using the relation %^ = F< x^ lrom the zc ~ 



< n k i > 



. k-l 



k\ 



{xf + sxD 



fcL ^ s!(fc _ s) , 

u>iC 2 (i) 
V 1 7TTWT ) < n li > + 

ad xvsw\ UiC 2 (i) 



+ {Xf - SX? W ) 

k-2 



+E 



k\ 



s\(k-s)\ 

s=0 v ; 



w» + Ci(i)k 
(_!)*-- 



< n s 2i > + 



+ 



+ d(i)k 
UiC 2 (i) 



LUi + Ci(i)k 



< > ] 



6 



< ni > 



fe-i 



A;! 



s=0 



(1 



s!(fc-s)! 
WjCi(i) 



(A7 d - « 2 st0 ) 



Wi + <7i(i)fc' " ' 22 



-{xi 

k-2 



5x sw 
k\ 



s=0 
[(1- 



S )! 

u)iC 2 (i) 



) < «5i > - 

Ui + Ci(i)k 
(_!)*-« 

+i 



< < > + 



+ C\{i)k J 



) < n s 2 p > + 



< n 



a+1 



> 



LUi + Ci(i)k 

The resulting equations for the zeroth moments couple 
to the higher moments by the interaction function F(i). 
These lower moments can be solved self-consistently. The 
resulting solution predetermines all the other moments, 
which completely describe the probability distribution. 
Each gene therefore couples to the other gene by the in- 
fluence of the self-consistently generated proteomic field. 
One could define the generating function and calculate 
the probabilics of having a given DNA binding state j 
for the i th gene when there are rii proteins of type i in 
the cell. In practice, it is easier to go back to the steady 
state master equation and solve directly for the proba- 
bility distributions than sum an infinite number of mo- 
ments. Rewriting the steady state master equation one 
gets: 



Xf d + 5Xr +Wi^^+n 



P 1 (n i =0) 



[(Xf + 6Xr)Pi(ni - 1) + 

+ + + 1) + UiFtim)] 

-[Pi(n i = l) 



Xf + SX? W + LU; 

LOiP 2 (ni = 0)] 



X ad _ 5X sw +LOl +n 

[{Xf-5Xt w )P2{^-l) + 

F(3 - i) 



P 2 (n t = 0) 



+(m + l)P 2 (n t + 
1 

Xf - 5X S ™ + lo, 
F(3 



"i 



1) +LJi 

[P 2 (n t : 
^Pi(n i = 0)] 



xl q 

1) + 



Pl{Ui 



These sets of equations give recursion relations for 
Pj(rii) which one can use to express Pj(n) as a func- 



tion of Pi(0) and P 2 (0)). The normalization condition 
En^o^iK) +P 2 {ni)) = l gives P,-(0) in term of con- 
stants and the result is the probability function Pj(n) 
as a series. The SCPF approximation reduces the two 
gene problem to a one gene problem parametrized by the 
moments of the second gene, which can be worked out 
independently, as we have already shown and are rep- 
resented by F(2), which is a constant in terms of this 
calculation. 

To see the effect of the stochastic nature of the system 
we compare the exact solutions of the self consistent field 
approximation equations to the results that would follow 
from deterministic kinetic rate equations for the number 
of proteins of each type and the fraction of on/off DNA 
binding states for each gene: 



Ci(l) 

Ci(2) 

n(l) 
n(2) 



X{ 



eg 



X 2 q 



X e 2 q +n2(l) 

Xf + 8X s 1 w {C l {l)-C 2 {l)) 
Xf + 5X s 2 w {d{2)-C 2 {2)) 



where n(i) is the number of proteins of type i present 
in the cell. The exact SCPF equations reduce to the 
deterministic kinetic equations in the limit of large u> and 
X ad for the case discussed above. The F(3 — i) term in 
the stochastic SCPF equations is replaced by the n 2 (3— i) 
term in the deterministic kinetic rate equations. For the 
toggle switch, where repressors bind as dimers it is easily 
shown that the interaction functional may be rewritten 
in the form: 

F{%) = (Xf) 2 +Xf + (5Xf) 2 + 

SXfidii) - C 2 {i)){\ + 2Xf) + 

-4uJi(5Xt w f 



^ 2 Ci(i)C 2 (i) 



uo t + Ci(i) 



=< n(i) > 2 



1 



LUi + Ci(i) 



-+ < n(i) > 



which in the large to limit reduces to F(i) =< n(i) > 2 
+ < n(i) >. So for large mean numbers of proteins 
present in the cell, which corresponds to large effective 
production rates X ad , < n(i) > of the order of hundreds 
is a small correction to < n(i) > 2 . We therefore repro- 
duce the deterministic kinetics result. 
As shown by Sasai and Wolynes (Sasai and Wolynes, 



7 




FIG. 2: Phase diagram obtained as an exact solution within 
the SCPF approximation for the single symmetric switch with 
X eq = 1 (a), 100 (b), 1000 (c). Contour lines mark values of 
AC. 




FIG. 3: Probability that genes are in the active state (a), 
the mean number of proteins of each type present in the cell 
(b) and the mean number of proteins present in the cell for 
gene i, when it is in the on state (c) as a function of X ad . 
Exact solutions of the SCPF approximation equations com- 
pared with deterministic kinetic rate equations solutions, for 
a single symmetric switch, X eq = 1000. 

2003) the difference in the probability that gene 1 is ac- 
tive and that gene 2 is active, AC = Ci(l) — Ci (2), plays 
the role of an order parameter. We can now consider a 
family of switches and discuss their stability, sensitivity 
of regions of bistability to control parameters and types 
of bifurcations. 



The Symmetric Toggle Switch 

For pedagogic purposes, we will start by analyzing the 
single symmetric toggle switch, such as discussed above 
in which repressors bind as dimers, with u>% = ll>2 = u>, 

X ad = X ad = X ad^ fijfsw = §X sw = §X sw &nd 

X 1 ^ 1 = X^ q — X eq , as it is the most intuitive and shows 
the most generic behavior. It is an academic example, 
as even individual genes in switches engineered in the 
laboratory mostly have different chemical parameters. 
Yet a lot can be learned from this simple system. 
The general mechanism of the phase transition. 



Figure [5] shows the phase diagrams for the system, 
| AC | as a function of reservoir protein number and the 
adiabaticity parameter for the exact SCPF equations 
for growing values of the parameter describing the 
tendency that proteins are unbound from the DNA, 
X eq . The deterministic kinetics and exact SCPF 
approximations give qualitatively similar results. The 
analogous deterministic kinetic phase diagrams agree 
with the SCPF solutions in the large u> and X ad limit, 
hence they become more similar with growing X eq , as 
the bifurcation occurs at larger effective production 
rates for larger X eq . For large fluctuations and a small 
unbinding rate, neither gene 1 nor gene 2 is favoured and 
the probability of a given gene to be on is determined 
solely by the effective production rate of the other gene 
and decreases in a quadratic manner as the number of 
repressor proteins grow (Fig. |3J. Since the switch is 
symmetric, the system has one stable state, AC = 0, 
where the probabilities of the genes to be on are equal. 
As the relative protein number fluctuations get smaller 
and the DNA unbinding rate grows, a proteomic cloud 
buffers the repressed gene, keeping it repressed. The 
symmetry of the system is broken and the solution 
bifurcates into two separate basins of attraction. For 
the stochastic SCPF equations the bifurcation takes 
place for larger effective production rates (larger X ), 
than for the deterministic equations, even in the large uj 
limit, which depicts their sensitivity to fluctuations. The 
critical number of reservoir proteins necessary for the bi- 
furcation of the solution to take place is the same in both 
approximations and is determined by < n > c — (X eq )2 
(Fig. |3J|. In the discussed example < n > c = 32 = 1000^, 
for X eq — 1000. For the deterministic kinetic switch the 
bifurcation takes place when C\{i) = — <n (3-i)>'i = i> 
due to the simple form of the interaction function equal 
to < n(3 - i) > 2 = (2X ad C 1 (3 - i)f . So d{i) = \ is 
equivalent to the < "^ e ^ > = 1. In a noisy system larger 
effective production rates are needed to achieve the criti- 
cal value of proteins. The interaction function in this case 
may be written as F(i) =< n(i) > 2 -^Q^m + < >> 



and 



;-W > 1, always. So at < n > c , 



> 1 



and the probability of the genes to be on is smaller 
than |, therefore C^ ff ' SCPF {i) < Cf ff ' kin {i). The 
mechanism of the bifurcation requires the two genes to 
be more likely to be unbound than bound for the phase 
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FIG. 4: Probability for each gene to be on (a) and mean 
number of proteins present in the cell (b) as a function of X ad , 
for different values of u = 0.0005, 0.005, 0.05, 0.5, 5, 50, 500, 
X eq = 100 for a symmetric switch. 

transition to take place. The curvature of the nullclines 
presented in Fig. [21 can be simply worked out to be of 
the form to = ^x^+^x^+b ~ 6, with Ct,£; constants 
determined by the specific value of Ci(l), Ci(2). 
Adiabaticity parameter dependence 
As the adiabaticity parameter decreases the area of 
phase space which corresponds to multiple solutions 
decreases (Figs |2 B}. For very small values of the 
adiabaticity parameter, there exists only one solution 
which corresponds to a state in which the two genes 
are off. The value of u> below which only one solution 
exists decreases with the tendency for proteins to be 
bound, but exists for all values of X eq . Therefore if the 
two genes have very high repressor binding affinities, 
the critical number of proteins necessary for the phase 
transition to take place cannot be formed, even for very 
high production rates. This region of parameter space 
where one solution is possible corresponds to a situation 
in which a buffering proteomic cloud may not form, due 
to a very fast destruction rate of proteins or a very small 
unbinding rate from the DNA. The critical number of 
proteins necessary for the bifurcation to occur grows 
with the tendency for proteins to be unbound from the 
DNA (X eq ), as the cloud buffering the genes needs to 
be bigger and exhibit smaller relative protein number 
fluctuations, which effectively decrease with the growth 
of the adiabaticity parameter. This is further discussed 
in terms of the probability distributions. Therefore a 
monostable solution exists at all values of the effective 
growth rate, X ad , for larger values of u> at large X eq than 
at smaller X eq values. The bifurcation point is a result 
of competition between the number of reservoir repressor 
proteins and the tendency for proteins to be unbound 



from the DNA. This is clear from the dependence of the 
number of proteins present in the cell at the bifurcation 
point on the relative values of X ad and X eq , but not the 
adiabaticity parameter lo, as can explicitly be seen from 
Fig. H 

Mean protein numbers 

The total number of proteins present in the cell, 
produced both in the on and off state, (Fig. 0J|, asymp- 
totically away from the bifurcation points is the same 
for the deterministic and stochastic approximations, and 
it is given by < n(i) >= 2X ad , when Ci(l) rj 1 the 
probability of the gene to be on is close to unity. The 
number of proteins of a given type present in the cell, 
when the gene that produces them is in the on state is 
always considerably smaller in the noisy system than the 
deterministic case (Fig. Efc)). Since the production rate 
in the off state was assumed zero, in the deterministic 
case no proteins of a given type are present in the cell if 
the gene is in the off state, unlike in the noisy system. 
Therefore the number of proteins in the deterministic 
system is nonzero only if the gene is on. But interaction 
of the DNA binding state with the proteins buffering it, 
results in a residual number of proteins present in the 
off state, for all values of lo. The region of bistability of 
the switch in parameter space grows as the binding rate 
increases with respect to the unbinding rate, stabilizing 
the DNA binding states. As the susceptibility of the 
system to fluctuations increases, the deterministic 
equations prove to be a poor approximation to describe 
the state of the system. 

Gene-buffering proteomic cloud interactions 

The stochastic nature of the system manifests itself 
also at the DNA level (Fig. 0). As the tendency for 
proteins to be unbound from the DNA grows, the area of 
parameter space, where multiple solutions are possible 
decreases, since a larger number of proteins is needed to 
reach a state in which two genes are more likely to be 
repressed (protein bound state), than at small X eq . For 
small unbinding rates or large binding rates, regardless 
of the ratio of the rate of unbinding of repressors from 
the DNA to protein degradation, bistability requires 
smaller numbers of proteins, which correspond to larger 
relative fluctuations, than for large X eq . Therefore 
a larger unbinding rate relative to the binding rate 
makes the system more susceptible to protein number 
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FIG. 5: Evolution of probability distributions for the proba- 
bility of the gene that will be active after the bifurcation to 
be on (a) and off (b ) as a function of the order parameter 
X ad . The bifurcation occurs at X ad = 44. 



FIG. 6: Evolution of probability distributions for the proba- 
bility of the gene that will be inactive after the bifurcation to 
be on (a) and off (b ) as a function of the order parameter 
X ad . The bifurcation occurs at X ad = 44. 




noise. Competition between X eq and < n(i) > results 
in X eq , for a given nullcline, being a parabolic function 
of X ad , for the dimer binding case, with coefficients 
determined by to and C\{i). This is easily generalized 
to higher order functions for higher order (p) oligomers, 
and results in p-order dependence. The switching 
region, by which we mean the region of parameter space 
between the bifurcation point and AC > 0.9 decreases 
as the binding and unbinding rates become comparable 
(X eq decreases). As discussed above, the probability 
of the genes to be on at the bifurcation point tends 
to i as the adiabaticicty parameter grows (Fig. 0}, 
therefore the probability to be on has to increase by a 
smaller AC to reach C\(i) = 1. Therefore the switching 
region decreases also as the unbinding rate from the 
DNA grows, since smaller effective production rates 
are needed to reach AC = 1, than for small u. Small 
values of co correspond to large fluctuations in the DNA 
binding state, as well as the protein number state and 
result in destabilizing the gene-buffering protein cloud 
interactions. Hence very large effective production rates 
are needed for AC > 0.9. Therefore the DNA unbinding 
rate must become considerably faster compared to 
proteins degradation rate for the switch to have two 
stable solutions in a large region of parameter space. 
The probability distributions 
A better understanding of the bifurcation can be gained 
from examining the probability distributions. Figures [S] 
and H3 show the evolution of the probability distributions 
of gene 1 and gene 2, respectively, to be on and off 
as functions of X ad . The peak of the distribution 
decreases and the width spreads out as the control 
parameter grows, until it reaches the bifurcation point 



at X ad — 44. Then the value of the probability function 
corresponding to the most probable number of proteins 
grows again. The spread of the functions grows as 
the effective production rate in the on state increases, 
however narrows with the increase of the adiabaticity 
parameter, as would be expected, since the DNA state 
fluctuations become smaller with u). The average number 
of proteins in the cell in the on state (AC > 0.9) does 
not show a dependence on uj. Yet as the unbinding 
rate from the DNA becomes very fast compared to the 
protein number fluctuations, the system switches often 
between the two states, hence a large number of proteins 
is present even in the off state. This results in a two 
peak - bimodal probability distribution (Fig. [5JH3). If 
the DNA unbinding rate is small, the protein number 
characteristics follow the DNA state having time to 
reach a steady state within each well, before the DNA 
binding site switches into the other state, so the number 
of proteins in the off state falls to zero (Fig. Q). If a; 
is large, random fluctuations in the DNA state do not 
change the effective state of the system, since a residual 
high mean protein number is present even in the off 
state. In such a case lower effective production rates 
than for small lj result in higher protein yields and what 
follows smaller switching regions. 

For small to one might expect Poisson distributions of 
proteins in each of the DNA states, since the unbinding 
rate from the DNA is smaller than the protein degra- 
dation rate, so the proteins may reach a steady state 
without the DNA state changing. Hence, effectively 
proteins would feel only one well and be subject to a 
birth death process. However this is not true. The 
difference between the exact solution and a solution 
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FIG. 7: Probability distributions for the gene to be in the on 
state (a) and off state (b) for a gene in the active state for 
different values of the adiabaticity parameter uo — 0.5, 10, 100. 
X eq = 100, X ad = SX SW = 100. 




FIG. 8: Comparison of probability distributions obtained by 
exactly solving the steady state equations in the SCPF ap- 
proximations with analogous Poissonian distributions. 

obtained within a Poissonian approximation to the state 
of the system is surprisingly large, owing to the skewed 
tails of these distributions. Figure |H1 compares these 
probability distributions with distributions for the same 
system if one assumes a Poissonian probability function. 
The distributions obtained as an exact solution within 
the SCPF approximation are clearly not symmetric, but 
exhibit long tails towards zero. Therefore, although the 
most probable values of the two types of distributions are 
similar, noise has a destructive impact on the system, 
resulting in a larger probability of having a smaller 
number of proteins in the cell than expected based on 
a Poissonian distribution, whose higher moments are 
equal to the mean. Therefore a larger production rate is 
needed for one of the states to be favoured as a result 
of noise than predicted from a symmetric probability 
distribution. The most probable number of proteins in 
the on state, if the unbinding from the DNA is slow, is 
zero, unlike predicted by Poissonian distributions. The 
influence of noise on protein number fluctuations brings 
the protein number means down, as can also be seen 
from Fig. [3] (c). Overall the spread of the probability 
distributions is large, and their characteristics for small 



FIG. 9: Nullclines for a symmetric switch when proteins bind 
as dimers for a symmetric switch when the effective base 
production rate </2/2fc ^ 0. Dependence on the adiabatic- 
ity parameter u) = 0.005, 0.05, 0.5, 5, 50 and the deterministic 
equations solution (a), <?2 = 5, X eq = 1000. Dependence on 
g 2 = 0.1, 0.5, 1.0, 5, 10, 10, lu = 0.5, X eq = 1000 (b). 
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FIG. 10: Probability of genes to be on (a) and mean number of 
proteins of a given type present in the cell (b) for a symmetric 
switch with an effective base production rate Sr = 5, uj — 0.5, 
x eq = 1000. 

values of the control parameters are different from 
those predicted by Poissonian distributions, let alone by 
deterministic kinetic equations, therefore the effects of 
stochasticity may not be neglected. 

The nonzero basal effective production rate 
case. 

The above analysis concerns a switch with a zero basal 
production rate, so proteins were not produced in the 




FIG. 11: Evolution of probability distributions for the prob- 
ability of the gene that will be active after the bifurcation to 
be on (a) and off (b ) as a function of the order parameter 
X ad for a system with a basal production rate ff = 5. The 
bifurcation occurs at X ad = 61. 
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off state. In a number of biological systems (Ptashne 
and Gann, 2002) a non-zero basal production rate exists 
and we now turn to consider the effect of this on a 
symmetric switch. Figure^ (b) shows the dependence of 
the bifurcation curves for different values of the effective 
basal production rate Values smaller than one, when 
the death rate is larger than the production rate, show 
that for the symmetric switch assuming the effective 
production rate to be zero in the off state is a reasonable 
approximation. If the on state has a positive input to 
the number of reservoir proteins present due to > 1 , 
the probability of the active gene to be on, even for very 
large on state effective production levels X ad is smaller 
than one. Hence the off state contributes considerably 
to the steady state number of proteins. The solution 
which corresponds to the more active of the two states 
may effectively be an off state, since it has C\{i) < ^, 
although the effective production rate in the on state 
in the bifurcated region of parameter space is much 
larger than in the off state (for example cyan line at 
|| = 20 in Fig. |^. As the effective basal production 
rate increases, a larger production rate in the on state 
than for small || > 1 is required to reach the critical 
number of proteins for the bifurcation to take place, 
which is given by < n(i) >= 2X ad d(i) -f(2d(i) - 1). 
For this reason even for the deterministic approximation 
at the bifurcation point, the two genes must be more 
probable to be off, as can also be seen for the exact SCPF 
solutions from the probability distributions fFigs lllU12|) . 
Figure El (a) shows the dependence of the bifurcation 
curves on the adiabaticity parameter, which tend to the 
deterministic case for large u>. A closer analysis of the 
> 1 case, since the 4j5f < 1 is analogous to the zero 
basal production rate case which was already discussed, 
show that mean properties of the system are in even 
better agreement with the deterministic solution than 
the c/2 = case (Fig. II Op . The system has a non-zero 
probability of being in the off state, with the probability 
distribution of the off gene having a long tail towards 
higher protein numbers Fig. El In the off state the 
effective production rate || is small and the noise input 
is small, relative to the large protein numbers present in 
the system. The small effect of stochasticity results in 
the observed similar mean characteristics. Yet the form 
of the probability distribution for the gene to be active 
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FIG. 12: Evolution of probability distributions for the prob- 
ability of the gene that will be inactive after the bifurcation 
to be on (a) and off (b ) as a function of the order parameter 
X ad for a system with a basal production rate 
bifurcation occurs at X ad = 61. 
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before the transition is especially broad, with far smaller 
probability than that of the inactive state (Figs 1111 112fl . 
These clearly show that the two genes are more probable 
to be in the off state before the bifurcation point. 
Therefore although the average observables are similar 
for the deterministic and SCPF stochastic solutions, the 
predicted distributions are unusual. 
Summary 

The symmetric switch is based on a competition between 
the accessibility of the repressor site and the number of 
repressor proteins present in the cell. The bifurcation 
is solely a result of the nonlincarity of the system and 
introducing noise simply affects the region in parameter 
space where given states occur. The protein number 
fluctuations have a destructive role in determining the 
stability of the bifurcated solution, however fast DNA 
unbinding rates can compensate for the destabilizing 
effect of protein number fluctuations. In this region 
the stochastic solution predicts similar means to the 
deterministic case, but the form of the probability 
distributions which depends on a large number of higher 
moments is non-trivial. It is a result of the interplay of 
the DNA binding and protein degradation kinetics. 



The Asymmetric Toggle Switch 

Most switches found in nature are not symmetric. For 
asymmetric switches, when proteins bind as dimers, the 
two genes interact, resulting in probabilities to be on, 
different from those imposed purely by the equilibrium 
between binding and unbinding. The steady state solu- 
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FIG. 13: Dependendce of probability of genes in an asymmet- 
ric switch to be on as a function of increasing parameters of 
one gene Xf d = 8Xf w (forward transition) for different values 
of X^\ 5 (a), 50 (b), 500 (c), keeping all other parameters 
fixed, X\ q = 1000, w x = w 2 = 0.5, Xf d = 5Xi w = 80, for 
deterministic and exact SCPF equations. 



FIG. 14: Dependendce of probability of genes in an asym- 
metric switch to be on as a function of decreasing parameters 
of one gene Xf d = 6X[ W (backward transition), for different 
values of X^ 9 : 5 (a), 50 (b), 500 (c), keeping all other param- 
eters fixed, X\ q = 1000, lo± = L02 = 0.5, X$ d = 6Xi w = 80, 
for deterministic kinetic rate and SCPF equations. 



MS (\ function of Ci(2) and: 



tion is a compromise between the tendency that repres- 
sors are unbound from the initially off gene {X^ q for the 
forward transistion, X^ for the backward in the follow- 
ing discussion) and the effective production rate of the 
initially on gene (X% d - forward, X^ d backward transi- 
tion) (at least for the deterministic case). This results 
in the characteristic S-curve bifurcation diagram, as pre- 
sented in, for example Fig. ^0 with possible forward 
and backward transtions, and what follows hysteresis. 
We refer to the transition which occurs with increasing 
Xf d as the forward transition and that with decreasing 
X® d as the backwards transition. Since Xf d is a well de- 
fined function of the probabilities that the genes are on, 
the simplicity of the deterministic equations allows for a 
completely analytic discussion of the asymmetric switch. 
The more complicated form of the exact SCPF equations 
makes this approach impossible. However the determin- 
istic rate solution offers valuable insight into the basic 
mechanism behind the transition. 
The general mechanism 
By combining the steady state equations of motion for 
the probabilities of the two genes to be on Eq. 2] and 
noting that with a zero basal production rate < n(i) >= 
2Xf Ci(i), one can derive the following form of the de- 
terministic bifurcation curves: 



Ci(2) 



1)< 



xf{c l {i)) = ^-{ 



2X% a 



1)5 (3) 



as a function of C\(l). The transistion points are de- 
termined as the extrema of these functions, which are 
functions solely of the scaled parameter X^/X^ 9 and 
are plotted on the bifurcation graphs. It is worth noticing 
that the bifurcation points C\{i) do not depend on the 
value of X^ 1 , the parameter describing the gene binding 
kinetics of the gene that is on initially. This is not true 
for the exact SCPF solution, which cannot be solved ana- 
litically, but the bifurcation curve has the more complex 
form: 

2 La(l) 1+oJi 
wi + <7i(l), 1 



(2) 



7 2Ci(l) 

where Ci(l) is a function of w 2 , X{ q , Ci(2) and X ad \ 
The bifurcation point is therefore determined by the 
protein [Xf ) and DNA (X^ 9 ) characteristics and their 
mutual interactions (u>i) of the two genes. The deter- 
ministic approximation therefore greatly simplifies the 
mathematical mechanism of the transition. This may 
lead to large errors when studying more complicated 
biologically relevant systems, where one considers asym- 
metric switches with non-zero basal production rates 
and proteins are produced in bursts. The case of the 
non-zero basal production rate within the deterministic 
approximation also cannot be solved analytically. 
The general picture behind the transition is seen from 
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FIG. 15: Evolution of the probability distribution for a for- 
ward transition as a function of X" d = 5X( W for X^ = 50 
with X\ q = 1000; Ul = u 2 = 0.5; Xf d = 6X$ W = 80 for an 
asymmetric switch. 




FIG. 16: Evolution of the probability distribution for a back- 
ward transition as a function of X* d = 5X™ for X^ q = 50 
with X\ q = 1000; ui = u 2 = 0.5; Xf = 8X% W = 80 for an 
asymmetric switch. 

the deterministic approach. The larger the tendency for 
proteins to be unbound from the DNA, the larger the 
effective production rate Xf d must be for the transition 
from one gene to be active to the other to be active 
to take place, since repressor proteins are less likely 
to bind to the on gene (i) at large X? 9 than at small 
X^ q . However, if one considers a noisy system, it is 
effectively harder for proteins to stay bound to the 
initially off gene due to the destabilizing effect of DNA 
binding noise (Figs El H4|l . For the stochastic system, 
apart from very low values of the adiabaticity parameter 
(u> < 0.1) (Fig. I19JI . there is a threshold number of 
reservoir proteins which will cause a rapid transition. 
If we start with a small effective production rate for 
one type of proteins and increase this rate, keeping the 
production rate of the other gene fixed at an initially 
higher value, the proteins produced by the gene with 
the initially smaller production rate, repress it gradually 
and ineffectively, until they reduce the probability of the 
gene to be on to one half, for the exact SCPF solution. 
The number of proteins present in the on state decreases 
much more rapidly with the change of Xf d , whether it 




FIG. 17: Mean number of proteins of each type present in 
the cell, according to exact solutions of the SCPF approx- 
imation for an asymmetric switch, with X^ 1 — 1000; uix = 
Lu 2 = 0.5; X 2 ad = <5XT = 80; and X% 9 = 5, 50, 500 during the 
forward (a) and backward (b) transition in the asymmetric 
switch. 

be increase for the forward transition or decrease for the 
backwards in the examples presented, than the number 
of proteins in the off state grows (Fig. I17JI . Hence the 
probability to be on of the initially active gene shows 
a larger sensitivity to the change of Xf d than the off 
state probability. This leads to a rapid transition of 
the previously active gene to an inactive state (Figs El 
llfijl . Such behavior is described by Ptashne (Ptashne, 
1992), (Ptashne and Gann, 2002) in the A phage switch, 
who points out its role as a "buffer against ordinary 
fluctuations in repressor concentration" . The observed 
system switches when the "repression probability" drops 
to 50%, as in the solutions of this model. Our analysis 
seconds Ptashne's hypothesis, since the deterministic 
system lacks this behavior, the transition is rapid and 
for certain values of parameters takes place when the 
probability of the initially on gene drops to 80% (Fig. 
I18fl . The buffering capabilities of the stochastic system 
are clearly seen in the long tails towards n = of the 
probability distibutions of the gene that is switching 
from the on to the off state (Fig. 115(1 . 
The effect of noise on the bifurcation mechanism 
The mean number of proteins at the transition point 
differs for the deterministic and exact SCPF so- 
lution (Fig. 117(1 . More repressors are needed to 
induce the transition in the deterministic approxima- 
tion than in the stochastic system, since due to the 
form of the interaction function for the exact case, 
F(i) =< n(i) > 2 zjf^+ < n(t) »< ni > 2 . A 
smaller number of proteins is therefore needed for the 
inactive gene to become competitive with the active 
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FIG. 18: Bifurcarion diagrams as a function of X^ d = 5Xf w 
for Xf q = 5, 50, 500 for exact solution of the SCPF and de- 
terministic kinetic equations for an asymmetric switch. Null- 
clines for Ci(l) with u> = 10 (a) and Ci(2) with lj = 0.5 (b) 
with X\ q = 1000;^! = lu 2 = 0.5; Xf = 5X% W = 80. 

gene. The mechanism of the transition is different from 
the symmetric gene case, where a critical number of 
proteins needs to be reached. The asymmetric switch is 
based on the competition between the probability that 
proteins of one kind will repress the opposing genes and 
the analogous probability for the other kind of proteins. 
The repression capability is governed by , which 
might be looked upon as the product of the probability 
of having a certain number of repressor proteins (3 — i) 
in the cell and the tendency for them to be bound to 
the opposing gene (i). In fact, the transition point 
in the deterministic case is purely a function of such 
ratios, £lj = f( x %n )■ In both the stochastic and 
deterministic cases, the transition points are set by the 
interaction function which regulates the on and off state 
probabilities of a given gene = §rjjy ■ Inclusion of 

noise in the system effectively increases the nonlinearity 
of the system, which results in the already discussed 
buffering capabilities of the system. Stochasticity alters 
the very simple competitive mechanism seen in the 
deterministic kinetics to allow for more levels of control 
of the stability of the state of the system against random 
fluctuations. 

Further comparison of solutions of the deterministic 
and stochastic equations leads to the same conclusions 
as for a symmetric switch. As the tendency for proteins 
to be unbound from the DNA grows, the difference in 
the critical number of reservoir proteins necessary for 
the transition to take place increases for both approx- 
imations. The critical number of proteins produced 
by a given gene necessary for the transition to take 
place for both genes is, in most cases (see u> dependence 
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FIG. 19: Bifurcation diagrams for an asymmetric switch, pre- 
senting Xt d = SXt w as a function of Ci(2) (top), and Ci(l) 
(bottom) for different values of the adiabaticity parameter: 
u>i=uj2 (a, e), 0J2, with u>i — 0.001 = const (b, d), u)i, with 
L02 = 0.001 = const (c, f). X\ q = W0,X^ q = 50,X$ d = 
5X$ W = 80. 

discussion), smaller for the exact solutions of the SCPF 

equations and the difference between the stochastic and 

deterministic result grows with both X^ q and decreases 

with oji (Fig. H8J. It has a value of 15 for X e 2 q = 500, 

cji = 0J2 — 0.5 and 2 for — 500, uji — UJ2 = 10. 

Consider the forward transition. The initially inactive 

gene is buffered by a cloud of repressor proteins. As one 

increases the effective production rate of the proteins 

produced by the inactive gene (Xf d ) : the number of 

proteins which are able to repress gene 2 grows slowly 

and linearly < n(i) >= 2Xf d Ci(l), where Ci(l) ~ const 

and forms a buffering proteomic cloud around it. In 

the results presented in the figures of this paper the 

tendency that proteins are unbound from gene 2, (X^), 

is smaller than A^ 9 , so gene 1 is able to produce enough 

repressors to form a stable buffering cloud around gene 

2 and turn it into the inactive state at quite modest 

values of X" d , If A^ 9 < X^ , gene 1 produces proteins 

less effectively, as the probability of it being repressed 

is larger than in the previous case, and larger values 

of Af d are needed to produce enough repressors to 

achieve a high effective probability of binding, y 1 ,, . An 

2 

example of how xf d ' crlt grows as A^ 9 — > X^ q , is seen by 
comparing the Xf - 33 for X{ q = 1000, A 2 e9 = 50 in 
Fig. [Hand Xf - 300 for X eq = 100, X eq = 50 (Fig. 
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Adiabaticity parameter dependence 

The interaction of the buffering proteomic cloud with 
the DNA can be altered when the rate of the DNA 
unbinding rate compared to the protein degradation 
rate is changed. For small u>i values the unbinding rate 
of repressors to the DNA is slower than the destruction 
of the produced proteins. Apart from very small u) 
values, as long as there is a critical number of repressor 
proteins in the buffering cloud, the off gene is repressed 
and it responds by turning on, only once the initially 
on gene is nearly totally repressed. Large adiabaticity 
parameters result in the efficient formation of the 
buffering proteomic cloud. For the initially off gene, a 
small DNA unbinding rate of the off gene, decreases the 
effectiveness of the buffering proteomic cloud around it, 
as the protein number state can reach a steady state 
before the DNA state does. The hindered DNA reaction 
to the protein number state effectively increases the 
tendency of repressor proteins to be unbound from 
the DNA, for a given X^ d . This in turn decreases the 
probability of the initially on gene to be on, leading to 
rapid, switching behavior as an be seen for gene 2 in 
the forward, or gene 1 in the backward transition for 
u> > 0.1 in Fig. E3(a). The initially on gene reacts to the 
interaction function of the initially off gene, for which 
F(i) ^>< n(i) > 2 + < n(i) > in the small id limit. 
Therefore the interaction function is effectively increased 
for C\{i) k, 0, leading to the enhanced buffering. The 
reaction of the initially off gene is unaltered, as for 
Ci(i) w 1 F(i) =< n(i) > 2 + < n(i) >~ const, if 
Ci(i) remains close to 1. However if u> is very small 
(black dash-dot curve in Fig. El (a)), the buffering 
proteomic cloud is not given a chance to form due to 
a very high degradation rate of proteins and gene 2 
is simply repressed in a gradual transition. If lj\ is 
extremely small and 0J2 large, the buffering proteomic 
cloud around gene 1 cannot form and the probability of 
it to be off in the forward transition decreases gradually. 
A buffering proteomic cloud exists around gene 2, hence 
the backward transition is reminiscent of the determin- 
istic result (Fig. El (b)). The most interesting case is 
shown in Fig. El (c), where a large io\ acts as a buffer 
against fluctuations in the number of proteins, which 
repress gene 1. For large production rates of repressors 



the probability of gene 2 to be on for the forward 
transition decreases faster than in the deterministic 
solution, however the buffering cloud repressing gene 1 
allows gene 2 to remain in the on state. A buffering 
proteomic cloud does not form around gene 2 and it 
remains on until the number of proteins produced by 
gene 1 grows considerably, as the effective production 
rate, A° d , is increased. The effective production rate of 
gene 1 must be very large to sustain a sufficient steady 
state number of proteins to repress gene 2 to the point 
that Ci(l) < 0.5, which leads to switching. For the 
backward transition the lack of a buffering proteomic 
cloud around gene 2 results in destabilizing gene 1 for 
larger Xf d effective production rates than for large u>2 
values. These examples show how certain combinations 
of values of adiabaticity parameters can lead to a system 
with a larger switching region than the deterministic 
model predicts. This property may be useful when 
engineering artificial switches. If one has a constraint on 
the production rates of the genes, one can use repressors 
with different binding affinities to achieve switching in 
the desired region of parameter space. 
In this simple system slow unbinding from the DNA 
can compensate for the destabilizing of the DNA state 
by protein number fluctuations. As the probability of 
the initially active gene to be on gradually decreases, 
the initially repressed gene becomes active only once 
the probability of the other gene to be on has fallen 
bellow a certain values a. The susceptibility of the 
system to protein number fluctuations may be estimated 
by the value of a. For small w, which is still able to 
sustain a buffering proteomic cloud, this values tends to 
0.5. The incapability of the system to form a buffering 
proteomic cloud is much stronger if both adiabaticity 
parameters are small, since the reaction of both genes to 
the change in the number of proteins is hindered (Fig. 
El (a)). DNA state fluctuations contribute to effectively 
faster protein number fluctuations, therefore the exact 
solution exhibits the very small u> characteristics, where 
a buffering proteomic cloud cannot form, for a slightly 
wider range of the adiabaticity parameter than one 
would expect with a Poissonian distribution (results 
not shown). Combining these observations a switch 
works most effectively if the change of the DNA state 
compared to the protein number fluctuations of one gene 
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FIG. 20: Bifurcarion diagrams as a function of Xf d = 5X{ W + 
2 f > with ^ = ^ = 5 (a) and = ^ = 0.5 (b) 

for X e 2 q = 5, 50, 500 for exact solution of the SCPF and de- 
terministic kinetic equations for an asymmetric switch. Null- 
clines for Ci(l) (a) and Ci(2) (b) with u)i — L02 = 0.5 with 
Xf = 1000; X$ d = 8Xi w = 80. 

is sufficiently smaller than that of the other gene, to 
allow for effective buffering. 
The nonzero basal production rate 

The asymmetric switch in which both genes have a 
nonzero basal effective production rate proves to be 
susceptible to noise. In Fig. |2S we show the dependence 
of Ci(l), with = ffigi = g and c ^ 2 ^ with 

g2 2 ^ = 92 S^ — 0.5 in the small u>i limit. The stochastic 
solutions converge to the deterministic solutions for 
large uj. If gene 2 is initially in the on state, the majority 
of proteins are produced with the high fixed rate in the 
on state, as g\ (2) >> 32(2). The repression of gene 2 
is in turn governed by the interaction function of gene 
1. If is small the number of proteins produced in 
the on and off states by gene 1 are comparable. As the 
number of proteins produced by gene 1 grows faster the 
larger g^. is, gene 2 gets repressed more effectively for 
smaller X^ d values. This results in a smaller number of 
repressors produced by gene 2 and the transition from 
gene 1 to be on to be off takes place for smaller Xf d - 
effective growth rate values, than for small 32- 
The deterministic solution is much more influenced by 
the production of proteins in the off state than the 
stochastic solution. In the exact SCPF solution slow 
DNA unbinding rates compared to protein degradation 
rates are another means of control of the stability of the 
DNA state against random protein number fluctuations. 
The state of the system is far less influenced by the exact 
protein numbers than in the deterministic solution. So 
until the probability of a gene to be on is larger than 
that to be off, the fraction of proteins produced with 



a smaller effective production rate in the off state is 
treated as a random fluctuation by the system. Once 
again the SCPF system demonstrates its susceptibility 
to protein number fluctuations. 

The influence of the off state protein production on 
the total repressor yield may also be seen in the fast 
decrease of C\{2) and increase of Ci(l) in the forward 
transition. If 52 is considerably large its effect can also 
be seen in the stochastic solution, hence even when gene 
1 is in the on state, it never reaches Ci(l) = 1, although 
gene 2 is totally repressed (Fig. |20l a and results not 
shown for gene 2). The magnitude of the probability 
of gene 1 to be on for very large effective production 
parameters strongly depends on the the tendencies of the 
proteins to be unbound from gene 1. As X^ increases 
the asymptotic X c { d limit of C\(l) becomes smaller, as 
it is effectively harder for repressors to stay bound to 
the DNA. The gene is more likely to be in the off state, 
which however manages to sustain the necessary number 
of proteins produced by gene 1 to repress gene 2. As 
g2 increases the region of bistability grows into areas 
of parameter space, in which the tendency of proteins 
to be unbound, X^ 1 , is larger than for small g^. For 
small values of X^ the number of repressors produced 
by gene 1 in the off state is sufficient to repress gene 2 
and one observes a smooth and slow transition in terms 
of Xf d . If 52 is considerably large the transition takes 
place for larger values of X^ d in the stochastic solution 
than in the deterministic solution, hence showing the 
large buffering region the interplay of DNA and protein 
number fluctuations provides. This also results in an 
effective similarity of the deterministic and stochastic 
solution. In regions of parameter space, in which the 
change of DNA state is rapid, the deterministic and 
stochastic solutions differ, apart from the large ui limit. 
Most experimentally observed proteins have very small 
basal production rates, which seconds our analysis, that 
it is functionally unfavourable for large basal production 
to occur. The dependence on other parameters is 
analogous to the case without a basal production rate. 
The region of bistability 
The backward transition, as already discussed, is analo- 
gous to the forward transition. In most cases, the regions 
of bistability (Fig. I18|l in parameter space are reduced 
in size by noise. When engineering artificial switches, 
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FIG. 21: Region of Ci(l) hysteresis for an asymmetric switch 
for the SCPF and deterministic approxiamtions as a function 
of u)i — u>2, with X^ q = 50 (a) and X^ q with uji — lu2 — 100 
(b). X{ q = 100, Xf = SXi w = 80, Xf = SX{ W . The 
same comparison for a switch with a basal production level 
H = 0.5 as a function of ui = u)i (c) and X^ (d). X^ q = 100, 
Xf = SXI W = 80. 

one may be interested in making sure the forward 
and backward transition takes place for considerably 
different production rates. We therefore consider how 
the region of bistability, defined as the difference in 
the critical effective production rate for the forward 
and backward transition, depends on the parameters 
of the model. For the deterministic case the region of 
bistability depends on the tendencies that proteins are 
unbound from the DNA in a quadratic manner, as can 
easily be seen from the bifurcation equations [3 El an d 
is demonstrated in Fig. The SCPF solution shows 
the same behavior. For large values of the adiabaticity 
parameter the size of the region of bistability is inde- 
pendent of w, as is the form of the bifurcation curve 
(Fig. ETTjl . The approach to this plateau is very rapid 
and is given by the ratio of polynomials. However, the 
size of the region of bistability for the u>i = ll>2 never 
reaches that of the deterministic solution, as even in the 
large u> limit the greater nonlinearity of the interaction 
function F(i) results in a more complex SCPF curve 
which does not reduce to deterministic solution, but 

xf{d(2)) -> micm - 1 ) x ?) h + ^ - l)5Z^iy * 
+ (2Xggi(2)) 2 ^_i_^ _ This effect ig tme 

for both curves, as the presented graphs show Ci(l) 



hysteresis and the chosen equations Ci(2). The same 
behavior is observed for the case with a nonzero basal 
production rate. The increase with X^ 9 is slightly slower 
in the 32 7^ case as the bifurcation curve is smaller by 
l^(Ci/(i)-C lin (i))-|/ng§|. 
Summary 

After the transition, the number of proteins produced 
by the now on gene, follows a linear dependence on 
X ad , similarly to the symmetric switch. The number of 
proteins in the cell is independent of the DNA dynamical 
characteristics, as those remain constant in that region of 
parameter space. The number of proteins of the off gene, 
rapidly falls before the transition takes place. Based on 
the bifurcation diagram of Fig. the phase transition 
is discontinuous, for a certain region of the parameter 
space, where switching may occur. That region may be 

roughly estimated by the parameters of the genes which 

x ad 9 x eq 

must be competitive, (-^3) ~ xf 7 ' This has a major 
implication for biological systems, such as the A phage, 
where many mechanisms are used to achieve balance 
between two genes. The first order phase transition, as 
opposed to the second order present in the symmetric 
system, is a results of the breaking of symmetry and is 
clearly seen in the evolution of probability distributions 
in phase space (Fig. I15|) . The gene that is on after the 
transition rapidly increases its probability of being on, 
whereas the off gene decreases with a rapid drop in the 
number of proteins it produces. 



The Case when Proteins bind as Monomers 

The equations presented above can easily be aug- 
mented to describe the binding of monomers or higher 
order oligomers by changing the form of the binding term 
to hifi^j, where p = 1 for monomers. The equations 
remain solvable for any value of p. 

Monomers do not make good repres- 
sors / activators 

The behavior of the system is quite different if we 
consider the case when proteins bind as monomers. For 
a symmetric switch there is no region of the parameter 
space, in which one observes switching. The SCPF equa- 
tions may be reduced to a single quadratic equation: 
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FIG. 22: Probability of genes in an asymmetric switch to be 
active when proteins bind as monomers, for different values 
of X e 2 q . X\ q = 1000,X? d = 5XI W = 80. 



2SX sw C 1 (if + (X eq +X ad -SX sw )C 1 {i)-X eq = (4) 

which has at most only one positive solution. Therefore 
the probability of one gene to be in the active state is 
always equal to that of the other to be in the active 
state and no switching is observed. The equation Q is 
independent of u, the adiabaticity parameter, therefore 
it is solely a consequence of the lack of nonlinearity 
in the binding of proteins and cannot be influenced 
by very slow DNA unbinding rates. By writing down 
deterministic equations we can also show that when 
proteins bind as monomers switching does not occur. A 
similar equation to Q), also independent of u>, holds for 
asymmetric switches. It also has one positive solution, 
therefore the parameters of the model predetermine 
the solution and each gene has a probability to be on 
determined by its kinetic rates. Since the rates are 
different for the two genes, the gene with the larger 
production rate will be in the active state, repressing the 
weaker gene (Fig. I22[) . In naturally occurring biological 
switches and those developed experimentally proteins 
bind as dimers, or higher order multimers (Ptashne, 
1992). We see cooperativity contributes to improving the 
efficiency of a switch. A switch controlled by monomers 
is shown to react ineffectively to changes in the repressor 
concentration, just as in the case of the asymmetric 
switch in our model discussed above. Monomers do 
not have the ability to stabilize a broken symmetry 
state, therefore the solution is fragile to kinetic rates 
and inefficient. Effectively monomers do not make good 
repressors/activators. Ptashne and Gann (Ptashne and 
Gann, 2002) explain the cooperativity process between 
two monomers by claiming that one monomer bound to 




FIG. 23: Probability distributions for the gene to be in the 
on state (a) and off state (b) for a gene in the active state for 
different values of the adiabaticity parameter u — 0.5, 5, 100. 
X eq = 1000, X ad = 8X SW = 50, when proteins bind as 
monomers to a symmetric switch. 

the DNA increases the "local concentration" of proteins 
around the binding site through weak protein-protein in- 
teraction, thus causing the second to bind cooperatively. 
Our model lacks spatial dependence, therefore shows 
this effect need not be thought of as due to changes 
in local concentration, but actually is required by the 
insufficient nonlinearity for monomers, which cannot 
produce bistability. 
Bimodal probability distribution 
Although the probabilities of the two genes to be on 
are equal for the whole region of parameter space and 
the mean number of both types of proteins in the cell 
is the same as in the deterministic case, the probability 
distributions are bimodal when the DNA unbinding rates 
are slower than the protein number fluctuations. The 
mechanism of this small u> behavior has already been 
discussed on the example of the symmetric switch when 
proteins bind as dimers. This is analogous to the case 
when DNA fluctuations induce a probability distribution 
with two peaks for the single gene with an external 
inducer (Cook et al., 1998). In fact the SCPF approxi- 
mation has reduced this two gene system to an effective 
one gene system with an external inducer. A bimodal 
distribution in the small u> case is also observed for 
the asymmetric switch, when proteins bind as monomers. 



The Case when Proteins bind as Higher Order 
Oligomers 

Switches in which effector proteins bind as higher 
order oligomers are omnipresent in nature and have been 
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FIG. 24: Phase diagram for the SCPF approximation for a 
single symmetric switch to which proteins bind as trimers (a) 
and tetramers (b), with X eq — 1000. Contour lines mark 
values of AC. 




FIG. 25: Mean number of proteins in the cell, for each type 
when proteins bind as trimers (a) and tetramers (b), ui — 
0.5, 10, X eq = 1000, symmetric switch. 

realized experimentally in artificial switches (McLure 
and Lee, 1998). We considered the binding of trimers 
(/ij(ri3_i) = hiTi^j) and tetramers (/ij(7i3_i) = /iinf_j) 
in symmetric switches. The equations of motion 
have the same form as before, but the interaction 
function F(i) accounts for the higher moments. For 
proteins binding as k th order oligomers it has the form 
F(i) = Ci(i) < n\(i) > +C 2 {i) < n\(i) >. As shown 
when discussing the dimer binding switch, the k th order 
moments have a simple form in the creation operator 
representation. 
The general mechanism 
From Fig. 1241 one notes that in order for the system to 
act as a bistable switch a considerably smaller number of 
reservoir proteins is needed than in the case of the dimer 
binding switch. As the multimericity number grows the 
area of bistability of the switch in parameter space grows. 
Since we assumed only one type of protein repressed 
a given gene, binding of higher order multimers is an 
effective model of cooperativity. Therefore we expect the 
system to have a larger region of bistability the higher 
the order of the binding multimer. The evolution of the 
system in parameter space when trimers bind is qualita- 



tively similar to the dimer binding scenario (Fig. I26fl . 
Fast DNA unbinding rates stabilize the system and the 
bifurcation takes place for smaller effective production 
rates, for large u> than for small u> (Fig. 125( 1. The critical 
number of proteins necessary for the bifurcation to take 
place is independent of the adiabaticity parameter and 
decreases with multimericity: < n > c — 32 for dimers 
binding, < n > c = 8 for trimers binding and < n > c = 4 
for tetramers binding. This along with the narrow 
probability distributions (Fig. I27JI . small u dependence 
when tetramers bind (Fig. I24JI . shows that one binding 
event determines the result, hence DNA binding rates 
do not play a role. Once there are < n > c proteins 
of a given type in the cell, a tetramer repressor will 
bind and stay bound. In the deterministic case the 
probability of the genes needs to fall to (p — l)/p, where 
p is the order of multimcrization of the repressor, for 
the bifurcation to take place. That along with the 
need for the number of repressors to be comparable 
with the tendency for proteins to be unbound from 
the DNA sets the critical number of proteins necessary 
for the bifurcation. Hence the bifurcation occurs when 
both genes are more probable to be on than off, for 
both tetramers and trimers. Therefore for the tetramer 
system a large buffering proteomic cloud is not needed 
to stabilize the DNA binding state of the switch and the 
characteristics of the system are practically independent 
of the adiabaticity parameter. 

Tetramer binding results in nearly deterministic 
characteristics 

In naturally occurring systems the production of the 
critical number of proteins is slowed down by relatively 
high multimerization rates and spatial dependence 
arising from the need of a large number of particles to 
diffuse together. These elements, which we neglect in our 
simple model constitute what might be called the cost of 
multimerization. This analysis also explains why most 
repressors and activators bind as dimers and tetramers, 
not trimers or pentamers. The effect of trimers binding is 
not different from that of dimers: a buffering proteomic 
cloud needs to be formed, the state of the system is quite 
influenced by noise, the switching region (region in X ad 
parameter space from the bifurcation point to AC > 0.9) 
is quite large. Yet in a real system there is an effective 
cost of trimerization: the energy of trimer formation and 
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all or none response to the protein distributions in the 
system. This formulation of the problem is naturally 
oversimplified, but it allows for general observations. 



FIG. 26: The evolution of the probability distribution for the 
probability of the gene that will be active (a) and inactive (b) 
after the bifurcation to be on as function of X ad for a switch 
when proteins bind as trimers, X eq = 1000, lj = 0.5. 




FIG. 27: The evolution of the probability distribution of the 
gene that will be active (a) and inactive (b) after the bifurca- 
tion to be on as a function of X ad for a switch when proteins 
bind as tetramers, X eq = 1000,w = 0.5. 



a need for the diffusion of particles. For tetramers the 
effect of stochasticity becomes negligible. Effectively one 
tetramer is sufficient for the bifurcation to take place. 
The binding of tetramer repressors may be thought of 
as a mechanism for increasing the deterministic nature 
of the switch. 

Binding of higher order oligomers as a com- 
petitive mechanism 

This analysis, although it neglects some important 
features, allows for a more quantitative formulation of 
cooperativity. Since most biological switches are asym- 
metric, cooperativity is also used as a means of making 
genes with smaller chemical rates more competitive. 
Tetramer binding seems to have a different role than that 
of lower order multimers. It may be used by genes which 
need to react to very small concentrations of proteins, 
for example they turn on degradation mechanisms when 
even a small number of toxic molecules is present. Or 
they may act as an extra mechanism stabilizing the 
existent state of a gene, as seems to be the case for the 
cl gene of the A phage. It seems tetramers are used 
as having either a stabilizing role or that of a drastic, 



The Case when Proteins are Produced in Bursts 



Many proteins in biological systems, for example the 
Cro protein in A phage are produced in bursts of N of 
the order of tens. We consider a symmetric switch, where 
proteins bind as dimers and are produced in bursts of N. 
The master equation in this case has the form: 



at 



dP 2 (n z ) 

at 



g 1 (i)[P 1 (n i -N)-P 1 (n i )} + 
+h[( ni + + 1) - mPxim 

-hinl_iPi(rii) + fiP 2 (rii) 
g2{i)[P2{n l -N)-P 2 {n l )] + 
+h[(ni + l)P 2 («-i + 1) - niP 2 (ni 
+h i n\_ i Px{ni) - fiP2(ni) 



for n > N. For n < N the equations have the form. 



dPijnj) 

at 



dP 2 (n l ) 

at 



-gi(i)Pi(m) + 

+h[(ni + i)Pi(n 4 + 1) - mPxim 

-hinl^Piin.,) + fiP 2 {ni) 
-9 2 (i)P 2 {ni) + 

+h[(ni + l)P 2 { ni + 1) - mPzirii 
+h i n\_ i Px(ni) - fiP 2 {rii) 



Following the same procedure as for the the single protein 
production case, we get the following equations of motion 
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for the first three moments 
dCi(i) 



at 
ac 2 (i 



at 

dC 1 ( t )<n 1 (i) > 

at 
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>]C 2 (i) 
> C 2 (i) 

>]Ci(») 

i) > C 2 (i) 
^N 2 ]C 2 {i) 
i) >]Ci(i) 
t) > C a (i) 



where = Ci(i) < n 2 (i) > +C 2 (£) < n 2 (i) > as be- 
fore. Writing out TV 2 = TV(TV— 1)+TV and subtracting the 
< rij(i) > equations from < n 2 (i) > we get the equations 
of motion for the previously defined creation operators a. 
Due to the form of F(i) for the dimer binding case only 
the first three moments are relevent. However generally 
this procedure can be carried out for higher moments, 
yielding an expression for the m th creation operator mo- 
ment in the steady state of the form: 



, 'lib - 

< a u . > 



{NXf + N5XD(l UiC2 £\.J < a?' 1 

y 1 A Lu l + mC 1 (i) J 1 





FIG. 28: Probability that gene i is on when proteins are pro- 
duced in bursts of TV — 10 (a) and TV = 100 (b), symmetric 
switch proteins bind as dimers, X eq = 100, u = 100. Com- 
parison of deterministic and stochastic solutions. 




(NX, 
TV™ -1 



ad 



nsxd 



uoi + mCi(z) 



'-{NXf -N8X? W (1 



- ra — 1 ^ 

< Qn > 



uJjC 2 (i) 
uji + mCi(i) 



FIG. 29: Mean number of proteins of each type present in the 
cell when proteins are produced in bursts of TV = 10 (a) and 
TV = 100, symmetric switch proteins bind as dimers, X eq = 
100, ui = 100. Comparison of deterministic and stochastic 
solutions. 



The general mechanism 

We discuss the effect of bursting phenomena on the ex- 
ample of a symmetric toggle switch when proteins bind 
as dimers, as that can offer the most insight, when com- 
pared to previous results. In this case switching takes 
> place for much smaller values of the effective production 
rate parameter X ad compared to when proteins are pro- 
duced separately. Therefore even in the large u> limit, 
noise resulting from large protein number fluctuations 
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To consider the binding of higher order oligomers when 
proteins are produced in bursts one simply accounts for 
the changed form of F(i) as discussed in the previous 
section. 



FIG. 30: Bifurcation curves as a function of X ad = SX SW , 
lo = 100 for different burst size values N = 1, 2, 5, 10, 50, 100, 
with X eq = 100 (a) and for proteins produced in bursts of 
TV = 100 (b) for different values of X eq = 1, 10, 100, 1000. 
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plays a role in denning the region of stability of the 
switch, as the criterion of large X ad is not reached. The 
number of proteins in the cell when the bifurcation oc- 
curs is determined by the tendency that proteins are un- 
bound from the DNA and does not change when pro- 
teins are produced in bursts. For the rates discussed in 
Fig. EH and Fig. ESI the critical mean number of pro- 
teins present in the cell at which the bifurcation occurs 
is n c = 10 = X eq = 100 5 . If proteins are produced in 
bursts of N = 10, as in the left hand figures, this value of 
n c is achieved when X ad > 1, that is proteins must get 
produced at a higher rate than they are destroyed to be 
able to sustain the steady state number of 10 proteins in 
the cell. In the figures on the right hand side of Fig. |2"5I 
and FigESl proteins are produced in bursts of N = 100. 
In this case even when the degradation rate is larger that 
the production rate, the critical steady state number of 
proteins necessary for the bifurcation to take place, can 
be reached and a bistable switch is possible. A bistable 
switch can exist even if the degradation rate exceeds the 
production rate for burst sizes present in biology. For 
X eq = 100, the order of the tendencies for proteins to be 
unbound from the DNA in the A phage, the value of N 
for which X" d < 1 is smaller than N = 20, the burst size 
for Cro proteins in the A phage. X ad at the critical point 
decreases as function of N (Fig. EOJ) and depends on the 
tendency that proteins are unbound from the DNA X eq 
(Fig. E01 (b)) and the adiabaticity parameter, uj (Fig. 

EUi. 

If proteins are produced individually the span of the non- 
adiabatic regime is clear from Fig. |2] It corresponds to 
w < 1. The bifurcation curves show small discrepancies 
for larger values of the adiabaticity parameter. However 
for larger burst sizes there is a continuous change in the 
form of the bifurcation curves with u>. All of the solutions 
differ substantially from the deterministic treatment, as 
shown in Fig. EH1 

The influence of the adiabaticity parameter on 
the bifurcation mechanism 

Contrary to the N = 1 case, the effective production 
rate at the bifurcation point X® d , grows smaller with 
the increase of the adiabaticity parameter, for consider- 
ably large burst sizes, as in the N — 100 example in 
Fig. El In this case each gene produces a large num- 
ber of repressors at a time. The bifurcation takes place 




FIG. 31: Bifurcation curves for proteins produced separately 
N = 1 (a), in bursts of TV = 10 (b) and N = 100 (c) as a 
function of X ad = 8X BW for different values of the adiabaticity 
parameter ui = 0.1, 1, 10, 100. 

in a region with X ad < 1, which corresponds to very 
small effective production rates, which denote very large 
death rates. Therefore in the region of parameter space 
before the bifurcation takes place both genes remain re- 
pressed (Ci(i) < 0.5) in the steady state, as opposed 
to the provisionally discussed situations, in which both 
genes had equal probabilities to be active {Ci{i) > 0.5). 
For large N bursts, the bifurcation takes place when one 
of the genes becomes unrepressed in the steady state. 
That is when the repressor cloud buffering the DNA be- 
comes destabilized, not when the cloud forms as in the 
smaller N examples. For large N bursts, if the rate of 
unbinding from the DNA is fast compared to the pro- 
tein degradation rate, larger effective production rates 
are needed for the buffering proteomic cloud to stabilize 
the DNA state, than for small u (Fig. E2 (c)). The 
larger X ad is, the more repressor molecules are present 
in the system, which corresponds to larger protein num- 
ber fluctuations, which are necessary for one of the genes 
to become unrepressed. For slower DNA unbinding rates, 
the buffering proteomic cloud is smaller, since the protein 
number reaches a steady state before the DNA state does. 
Therefore the buffering proteomic cloud is destabilized at 
smaller values of X ad . Hence, in the case of small uj the 
unrepressing bifurcation takes place for smaller effective 
production rates than for large uj. However if the un- 
binding rate from the DNA is very small, u> < 0.01, X^ d 
as a function of the adiabaticity parameter grows again, 
as this corresponds to effectively large death rates, which 
need very high production rates to sustain a proteomic 
cloud buffering the DNA. If the effective production rate 
is too small in this case, the steady state number of pro- 
teins is too small to form the buffering proteomic cloud, 
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FIG. 32: Bifurcation curves for proteins produced in bursts 
of N = 30, 55, 56, 57, 100 as a function of X ad = SX SW for 
different values of the adiabaticity parameter to = 0.1, 100. 

although the burst size is enormous. In the very small 
u> limit the bifurcation cloud needs to be formed for the 
bifurcation to be possible, as in the mechanism present 
in the small N case. The value of X ad at the bifurca- 
tion point in both the large and small to limit is strongly 
governed by protein and DNA binding state fluctuations 
in the system. For this reason the deterministic solu- 
tion fails. It assumes the incorrect mechanism, in which 
the bifurcation is a result of repressing one of the genes. 
This can happen if the death rate of proteins is slow 
enough to allow for the existence of < n(i) c > repressor 
molecules in the system at very small production rates 
(Ci(l) h4// < fcm = 0.5) (Fig. HHJl. One can see that the or- 
der of taking the adiabatic limits in the steady state for 
proteins produced in large bursts is subtle and depends 
strongly on the parameters of the system, as the bifurca- 
tion is governed mainly by relative protein and DNA fluc- 
tuations, both of which are very large. Furthermore, the 
deterministic solution is closer to the small u limit, which 
corresponds to slow DNA unbinding rates compared to 
protein number fluctuations. Deterministic results may 
therefore be misleading in the bursting situation, even 
for large uj. 

Figure 02 shows explicitly how the steady state comes 
about as a result of different mechanisms depending on 
the burst number TV and how the order of reaching the 
steady state by the protein and DNA binding site dynam- 
ics changes depending on u>. The other parameters were 
chosen so the bifurcation would take place for X ad < 1 
for all examples. For small burst sizes, slower DNA un- 
binding rates require larger effective production rates to 
reach the steady state number of proteins necessary to 
form the buffering proteomic cloud than for large N . For 
larger burst sizes, faster DNA unbinding rates destabilize 




FIG. 33: Probability for genes to be on in the SCPF ap- 
proximation as a function of X ad = 5X SW , with X eq = 100, 
u = 100 for N = 1, 2, 5, 10, 20, 50, 100 (a) and the mean num- 
ber of proteins present in the cell (b).) 

the buffering cloud of proteins for smaller effective pro- 
duction rates than in the small N case. The inset shows 
the transition region between the two possible mecha- 
nisms. 

Consequences of bifurcation at smaller X ad val- 
ues 

The divergence from the deterministic solution at the 
bifurcation point increases with the burst size, as is ex- 
pected due to the enormous noise effect due to large N, 
on a system with a constant and independent of the burst 
size number of proteins at the bifurcation point (Fig. l3"3"|) . 
As already noted the number of proteins in a cell, is 
in the range of tens to hundreds, even if they are pro- 
duced in bursts. Figure 1331 shows that this number is 
reached for smaller effective production rates for larger 
burst sizes than for small N values. Therefore systems 
where proteins are produced in bursts display smaller 
values of X ad and are more susceptible to noise if the 
number of proteins in the cell is to be of the order which 
is observed experimentally. Furthermore the noisy burst 
systems even for very large values of X ad do not con- 
verge as closely to the deterministic solution as they do 
for the single protein production example. This can be 
seen from the form of the steady state moment equations. 
The interaction function F(i) for the N = 1 case in the 
limit of large u> and X ad converges to F(i) — >< n(i) > 
+ < n(i) > 2 whereas the deterministic solution corre- 
sponds to F(i) =< n(i) > 2 . Therefore for large mean 
values of proteins the two are equal. However in the case 
when N > 1, F(i) ^< n(i) > (1 + E f 1 )+ < n(i) > 2 , 
which requires N << 2 < n(i) > for the effect of burst- 
ing to be negligible at very large N. The values of the 
effective production rate that correspond to values of the 
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proteins seen experimentally seem to be small. There- 
fore we can say that effectively the role of bursting is to 
enable for the existence of a bistable solution at lower 
effective production rates, which determines a region of 
parameter space which has been previously unstudied. 
In this region one cannot make the adiabatic assump- 
tion that the change in the DNA state can be integrated 
out due to a separation of timescales. That assumption 
leads to erroneous results, predicting a region of bistabil- 
ity where explicit treatment of both timescales suggests 
monostability. Furthermore, for very large TV, the region 
of bistability decreases with the adiabaticity parameter, 
making the disagreement of the stochastic solutions with 
those of the deterministic rate equations larger. The adi- 
abatic approximation and the full solutions converge only 
in the regime of large u> and X ad , the second of which 
is never fulfilled at the bifurcation point or for biological 
concentration for systems in which proteins are produced 
in large bursts. 

Dependence on the DNA Binding Coefficient 

Just as increasing the burst size, decreasing the tendency 
for proteins to not be bound to the DNA results in a 
different switching mechanism. The probability of the 
genes to be on falls to far smaller values than the 0.5 
of the N = 1 case. If the burst size is large both genes 
have a very low probability of being on before the critical 
number of proteins necessary for bifurcation is achieved. 
The same effect is observed if proteins are more likely 
to bind to the DNA (small X eq ) (Fig. I3UI (b)). When 
the genes are more probable to bind a repressor and suc- 
cessful unbinding events are rare, earlier bifurcations in 
terms of X ad result. As X eq increases, the probability of 
the genes to be on at the bifurcation point decreases as 
repressors have a higher tendency of unbinding. 
For very high values of the adiabaticity parameter, corre- 
sponding to high unbinding rates form the DNA binding 
site, the stable solution which corresponds to the off state 
and the unstable state merge and the system is monos- 
table again, with only the on state present. This limit is 
also reached by keeping X ad fixed but taking the burst 
size TV — > oo. 
Probability distributions 
In the case of the rates used in Fig. [M] and Fig. [351 
n c = 32 is the same as for TV = I, but we note a tenfold 
decrease in X^ d compared to when proteins are produced 




FIG. 34: The evolution of the probability distribution of the 
gene that is on after the bifurcation, to be on (a) and off (b) 
as a function of X ad for a switch when proteins are produced 
in bursts of TV = 10, X eq = 1000, cu = 100. Bifurcation point 
at X ad = 8X SW = 35. 




FIG. 35: The evolution of the probability distribution of the 
gene that is off after the bifurcation, to be on (a) and off (b) 
as a function of X ad for a switch when proteins are produced 
in bursts of TV = 10, X eq = 1000,oj = 100. Bifurcation point 

at X ad = 5X SW = 35. 

separately. When proteins are produced in bursts, the 
probability distributions have tails towards larger n, as 
opposed to the distributions for individual protein pro- 
duction. The mean number of proteins in the system 
for given states of the switch is similar to that of the 
N = 1 case, however the distributions with bursts are 
much broader, as could be expected. In this case even 
very fast unbinding rates from the DNA cannot correct 
for the enormous protein number fluctuations and one 
must explicitly keep track of the change of the DNA bind- 
ing state. A system in which proteins are produced in 
bursts is very noisy, especially compared to the nearly 
deterministic case of proteins binding as tetramers. 

Nonzero basal effective production rate 
If there is a nonzero basal production rate the difference 
between the deterministic and stochastic solutions is also 
qualitative even for relatively small burst sizes. In this 
case proteins are also produced in the off state, so there 
the number of repressors produced by the off gene after 
the bifurcation is nonzero, but equal to the burst size TV, 




FIG. 36: Probability that gene i is on when proteins are pro- 
duced in bursts of TV = 10 with a basal effective production 
rate g 2 /2k = 0.5 (a) and TV = 100, with a basal effective pro- 
duction rate g 2 /2k — 0.05 (b), symmetric switch proteins bind 
as dimers, X eq = 100, uo — 100. Comparison of deterministic 
and stochastic solutions. 




FIG. 37: Mean number of proteins produced by each gene 
when proteins are produced in bursts of TV = 10 (a) with a 
basal effective production rate g 2 /2k = 0.5 and TV = 100, 
with a basal effective production rate g 2 /2k = 0.05, sym- 
metric switch proteins bind as dimers, X eq — 100, u = 100. 
Comparison of deterministic and stochastic solutions. 

since < n(i) >= N(X ad + 5X SW (2C 1 (i) - 1)) -^W^ 
N2%r. This number is equal for both the stochastic and 
deterministic solutions and is equal to 10 in the examples 
presented in Fig. 03 So production in bursts maintains 
a high level of repressor proteins, even for very small 
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FIG. 38: The evolution of the probability distribution of the 
gene that is on after the bifurcation, to be on (a) and off (b) 
as a function of X ad for a switch when proteins are produced 
in bursts of TV = 10 with a basal effective production rate 
g 2 /2k = 0.5, X eq = 100, to = 100. bifurcation point at 
X ad = SX sw + 2 g 2 /2k = 8. 
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FIG. 39: The evolution of the probability distribution of the 
gene that is on after the bifurcation, to be on (a) and off (b) 
as a function of X ad for a switch when proteins are produced 
in bursts of TV = 10 with a basal effective production rate 
g 2 /2k = 0.5, X eq = 100, lu = 0.5. Bifurcation point at X ad = 
SX BW +2g 2 /2k = 6. 

% values if the burst size is large. When using exper- 
imental data one must be very careful to consider the 
burst size when assuming the basal production level is 
zero. Furthermore, the value of the interaction func- 
tion of the gene in the off state (Ci(i) ~ 0) for the 
stochastic case is much larger than for the deterministic 
case, due to the multiplication of < n(i) > 2 which gives 
F(i) -K n(i) > 2 (1 + 5^) + N%, for large u, the effect 
of which is shown in Fig. 1361 The number of repressor 
proteins produced by the off gene decreases as g2 — > 0, as 
expected and the probability of the on gene to be active 
tends to one, as is shown in Fig. 0U](a). The dependence 
of the effective production rate at which the bifurcation 
occurs on the adiabaticity parameter is analogous to that 
of 52 = case. The very small ui cases are shown explic- 
itly in Fig. The probability distributions for the gene 
which is active after the bifurcation in the on and off state 
are presented in Fig. 1381 for large unbinding rates from 
the DNA, and Fig. 023 for small unbinding rates from the 
DNA. They exhibit maxima around 2X ad for the on state 
and 2|r for the off state and display behavior analogous 
to that of proteins produced separataly, apart from the 
different curvature of the slopes for n < TV and n > TV. 
For small uj values the protein numbers reach a steady 
state before the DNA states, hence we observe bimodal 
probability distributions. The mechanism of competition 
in this noisy burst system is different than in the single 
protein production case. If the gene is in the on state, 
probability states with higher n values are strongly oc- 
cupied and there is hardly any probability flux into the 
lower n states. In the off state however, a flux pushes 
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FIG. 40: Bifurcation curves as a function of X ad = M™, 
X eq = 100, lj = 100, N = 100, for different basal ef- 
fective production rate g 2 /2k = 0.0005,0.005,0.05,0.5,1,0 
(a). Comparison of lo dependence with deterministic so- 
lution. N = 20, g 2 /2k = 0.5, X eq = 100, uj = 
0.0005, 0.005, 0.05, 0.5, 5, 50 (b). 

the system into the lower n states, essentially trapping it 
there, hence the difference in the slopes, as can be seen 
in Fig. |2H This is also true for the g2 — system when 
proteins are produced in bursts. 



Limitations of the SCPF Treatment 

The examples presented above cover a large class of 
two gene switches, all of which are exactly solvable within 
the SCPF approximation. An exact solution may be ob- 
tained within this approximation for systems of genetic 
networks and switching cascades. However the SCPF ap- 
proximation does not allow for an exact analytical solu- 
tion of all systems. If wc try to model one of the simplest 
natural systems where regulation is achieved by means of 
a switch, that is the A switch, we encounter a problem. 
The genes in the A switch, apart from having a toggle like 
regulation, also exhibit auto-regulation, that is cl pro- 
teins can bind to OR3, repressing the cl gene, and the 
Cro proteins can bind to OR1 or OR2, enabling the RNA 
polymarase from transcribing the Cro gene (Ptashne, 
1992), (Ptashne and Gann, 2002). If we expand the mas- 
ter equation to account for self-regulation we add a /i^nf 
binding term to the -Pj(nj) equations. Therefore the k th 
moment equation will display a dependence on the k+p th 
moment and the set of equation will not exhibit closure. 
One can find the probability distribution for a single self- 
regulating single gene. However if we consider as system 
like the A phage, where self regulation is also combined 
with regulation by another gene, the problem is no longer 



solvable exactly and demands a cutoff of the hierarchy or 
other approximations. We can nevertheless treat these 
systems using the variational method, as proposed by 
Sasai and Wolynes (Sasai and Wolynes, 2003). The fact 
that self-regulation renders the system incompletely solv- 
able within the SCPF approximation, is not surprising, 
since it corresponds to the exact solution for such a sys- 
tem. Gene i is influenced only by the number of proteins 
it produces. It is independent of the state of the other 
gene. Therefore, as one would expect the full solution 
should depend on all moments of the distribution of gene 
i. However for systems such as the A phage, we can treat 
all inter gene regulation effects exactly and truncate the 
self-regulation equation at the highest order of the inter 
gene interaction, which would be six, corresponding to, 
for example, 3 cl proteins binding to the 3 operator sites. 

Conclusions 

The self-consistent proteomic field approximation for 
stochastic switches reproduces many intuitive notions 
about their behavior. It proves to be a a very power- 
ful tool that allows for the consideration, of all but one, 
of the basic building blocks of more general switches and 
networks. A switch with a self-repressing/activating gene 
cannot be solved exactly within the SCPF approxima- 
tion, as in this case the approximation is equivalent to 
the full solution. Therefore the probability distribution 
is determined by an infinite number of moments. The 
probability distributions obtained for the systems con- 
sidered in this paper are not symmetric and exhibit long 
tails. This anticipates problems for using the variational 
principle for finding probability distributions when one 
accounts for correlations between the two states. The 
possibility to expand this method to consider networks 
and cascades will allow for are more realistic treatment 
of complex systems with emergent behavior at low com- 
putational costs. 

One can account for the mRNA step in the system by 
a adding a deterministic step which using a determinis- 
tic kinetic rate equation translates the number of mRNA 
molecules into proteins produced in bursts. This is a 
valid procedure, as as separately shown by (Thattai and 
van Oudenaarden, 2001) and (Swain et al., 2002), tran- 
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scription noise is just amplified in the translation process. 
Therefore treating the mRNA step deterministically sim- 
ply introduces another constant into the discussed case 
of proteins produced in bursts. Therefore the presented 
treatment of proteins produced in bursts with a modified 
effective production rate is a simple model of including 
mRNA in the system. Of course, the effect of mRNA is 
much more complicated, as it also introduces, for exam- 
ple time delay, between binding and production. This 
model in the present state neglects these effects. 
Our analysis of a large class of switches, shows how par- 
ticular elements contribute to the emergent behavior of 
functioning switches. Comparison of the stochastic and 
deterministic treatments of a single gene switch shows 
convergence in the region of fast rates of unbinding from 
the DNA compared to protein number fluctuations and 
large effective production rates. For symmetric switches 
when proteins are produced separately the two solu- 
tions converge after the bifurcation, but often differ when 
defining the region of parameter space, where the bifur- 
cation occurs. The agreement between the deterministic 
and stochastic solutions, is especially good for symmet- 
ric switches, with N = 1 and a non-zero basal production 
rate. However even though the mean repressor protein 
levels in the cell are similar in both approximations, the 
probability distributions arc broad and far from Poisso- 
nian, i.e. they are not completely characterized by these 
means. If the adiabaticity parameter is small (u) < 1) 
the protein number state reach a steady state before the 
DNA binding state and we observe a bimodal probabil- 
ity distribution. For the symmetric switch noise has a 
destructive effect on the region of bistability. Increas- 
ing the adiabaticity parameter facilitates the formation 
of a buffering proteomic cloud around a gene, which leads 
to repression at lower effective production rates than for 
small uj. 

As was already mentioned, the symmetric switch is hard 
to design and build experimentally. The asymmetric 
switch, which is the experimental toy system, is much 
more susceptible to noise than the symmetric switch 
and stochasticity has not only the destructive effect on 
the region of stability one might expect, but also intro- 
duces new phenomena and can be utilized to increase the 
bistable region. This is of fundamental importance, since 
experimentally one deals with asymmetric switches and 



these offer greater possibilities in artificially engineering 
new systems. As can also be learned from the asymmetric 
switch as well as from the analysis of binding of differ- 
ent oligomers, the region of bistability of a switch grows 
with increasing the interaction function. When creating 
artificial switches, one may argue a large region of bista- 
bility may be desired, so the switch reacts by the forward 
or backward transition to very specific concentrations or 
production levels of a protein. If the experimental setup 
constrains the protein production rates, this can also be 
achieved by modifying the adiabaticity parameters of the 
system, which ensures the transition remains rapid and 
effective. Asymmetric switches, exhibit first order phase 
transitions. This size of the region of phase space, in 
which the forward and backward transitions occur grows 
with the tendency that proteins are unbound from the 
DNA of both genes. Large adiabaticity parameters sta- 
bilize the buffering proteomic cloud around the repressed 
gene and lead to the formation of an effectively repressing 
cloud for smaller numbers of repressors, in the forward 
transition, than for small ui, for the active gene. 
Experimental data available at this point (Darling et 
al., 2000), suggest biological switches function in regions 
of high adiabaticity parameters from the deterministic 
point of view. Nevertheless, even for large values of 
adiabaticity parameters one must account for the DNA 
binding site fluctuations explicitly when proteins are pro- 
duced in bursts. The deterministic solutions give quali- 
tatively wrong results in biologically relevant areas of pa- 
rameter space. The stochastic solutions for large burst 
sizes suggest that the bifurcation of the solution is a re- 
sult of destabilizing of the repressor cloud buffering the 
DNA, not formation of the cloud as for smaller burst sys- 
tems. The probability distribution therefore exhibit tails 
towards large n values, not as in the small N case to- 
wards small n values. The deterministic kinetics remains 
unchanged for large burst sized, unlike the stochastic 
kinetics, hence presenting results derived from a wrong 
mechanism. The definition of the adiabatic limit, when 
proteins are produced in bursts is not clear as in the 
N = 1 case, when it corresponds simply to w < 1. This 
ambiguity does not allow one to integrate out the degrees 
of freedom corresponding to the change in DNA binding 
site occupation. Such an approximation leads one to er- 
roneously identify the regions of bistability. The switch 
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with a nonzero basal production rate when proteins are 
produced in bursts results in probabilities to be on and 
mean numbers of proteins in the cell very different from 
those of the deterministic solution, even for small effec- 
tive basal production rates. If proteins are produced in 
bursts assuming that a small effective basal production 
rate may be approximated by a zero rate may be mis- 
leading. Binding of proteins produced in bursts results 
in a bifurcation transition for smaller values of the effec- 
tive production rate. It is also a mechanism for making 
two genes in an asymmetric switch more competitive. 
Binding of higher order oligomers leads to results closer 
to those of deterministic treatments, with narrower prob- 
ability distributions. This can be experimentally used 
to stabilize DNA binding states. In this simple model 
tetramers seem to be the most optimum binders,. The 
close to deterministic all or nothing switching they offer 
may be worth the effective cost of the energy of multi- 
merization and diffusion of particles. Binding of higher 
order oligomers may be viewed as a simple model of coop- 
erativity, which increases the competitiveness of genes in 
an asymmetric switch. Within the SCPF approximation 
monomers do not make good switches due to lack of non- 
linearity in protein concentration. They do not exhibit 
a region of bistability. This model neglects any struc- 
tural DNA-protein interactions and spatial dependence. 
Hence this conslusion is simply a result of the lack of 
cooperativity in the system. For small adiabaticity pa- 
rameters, they do however exhibit bimodal probability 
distributions, unlike in the large ui limit. 
The thorough investigation of different components of 
gene regulatory networks using the self-consistent pro- 
teomic field approximation provides a tool kit for engi- 
neering new switches and networks. Based on our anal- 
ysis, if one would want to build a strong component of 
a switch out of a gene with relatively small chemical pa- 
rameters, one could use components that utilize binding 
of tetramers and that produce proteins in bursts. This is 
what the Cro gene in the A switch uses. 
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